Graphene¶

Graphene consists of carbon atoms in a 2d hexagonal configuration. Like WSe$_2$, graphene is a popular material for twisting and can be considered the cradle of twistronics since Macdonald’s seminal publication. [BM11] Due to its simplicity and similarity to our research interest, it is a good example to showcase the application of Quantum Espresso (QE) [GAB+17] and Wannier90 (W90) [PVA+20] to a 2d material.

Basic theory¶

A possible choice of lattice vectors for graphene is:

(21)¶\[\begin{split} \vec a_1=\begin{pmatrix}1\\ 0\end{pmatrix}\text{ and } \vec a_2=\frac12\begin{pmatrix}-1\\ \sqrt{3}\end{pmatrix}, \end{split}\]

where length is in units of the lattice spacing $a$ between two carbon atoms, which is approximately 1.42 Å. [CDAnjouG+11]

../../_images/logo2.png — Fig. 13 Lattice structure of graphene.¶

In the case of sublattice symmetry the Hamiltonian writes:

(22)¶\[\begin{split} H(\mathbf k) = \begin{pmatrix} 0&h_\mathbf k\\ h_\mathbf k^\dagger&0 \end{pmatrix}, \end{split}\]

with:

(23)¶\[ h_\mathbf k = 1+\text e^{-i\mathbf k\cdot\vec a_1}+\text e^{i\mathbf k\cdot\vec a_2}, \]

where energy is measured in units of the hopping parameter $t\approx 2.7$ eV between two carbon atoms [CDAnjouG+11].

In Fig. 14 we plot the band structure for graphene for this tight-binding Hamiltonian. At the K point the dispersionbecomes linear and forms a Dirac point.

Quantum Espresso¶

Now we wish to recreate this bandstructure and Hamiltonian with a dft calculation. The first step here is a self-consistent field calculation which fixes the electron density for future band structure calculations.

SCF/NSCF/Bands input¶

scf &CONTROL input

&CONTROL
  prefix = 'graphene'
  outdir = './out'
  verbosity = 'high'
  pseudo_dir = '../pseudo/'
  calculation = 'nscf'
/

Description

The &CONTROL card contains parameters controlling IO

QE finds the required pseudopotentials in the pseudo_dir and stores them and all datafiles in the outdir. The prefix functions as a label for these files so that later calculations can find these specific files and reuse them.
verbosity sets the amount of information to be written to the output file.
calculation type, in this case scf.

scf &SYSTEM input

&SYSTEM
  assume_isolated = '2D'
  ibrav = 4
  nat = 2
  ntyp = 1
  occupations = 'smearing'
  degauss = 0.015
  ecutwfc = 45
  nbnd = 30
  a = 2.45951214
  c = 32
/

Description

The &SYSTEM card sets parameters controlling lattice properties

QE contains a number of preset lattices, with ibrav 4 we choose a triangular lattice with lattice vectors:
1. $v_1 = a\begin{pmatrix}1&0&0\end{pmatrix}$
2. $v_2 = a\begin{pmatrix}-\frac12&\frac12\sqrt3&0\end{pmatrix}$
3. $v_3=a\begin{pmatrix}0&0&\frac ca\end{pmatrix}$
Here a and c controll the size of the lattice. For 2d lattices, like graphene, we use assume_isolated and a large value for c to make sure there is no periodicity in the z direction. Finally nat and ntype inform QE of the total number and the number of unique atoms in our unit cell.
In the case of metals DFT can be unstable while integrating the Fermi surface. We prevent problems by setting occupations to smearing which smears the bands with a value set by degauss.
ecutwfc is a cut-off for the plane wave expansion of the wavefunctions. We can find appropriate cut-offs for each atom type at the materialscloud. [PMC+18]
nbnb is the number of bands we wish QE to calculate. For a complete picture of the conduction bands it is important to set this parameter manually.

scf &ELECTRONS input

&ELECTRONS
/

Description

The &ELECTRONS card must be included for a scf calculation even though we use only the default parameters.

scf ATOMIC inputs

ATOMIC_SPECIES
  C    12.0107   'C.pbe-n-kjpaw_psl.1.0.0.UPF'
ATOMIC_POSITIONS angstrom
  C    0.0000000   0.0000000   0.0000000
  C    0.0000000   1.4200000   0.0000000

Description

The ATOMIC_SPECIES card is a list of all atom types with their weight and pseudopotential files, non-relativistic pseudopotentials can also be found on the materialscloud. [PMC+18]
The ATOMIC_POSITIONS card sets the unit of length and the positions of the atoms in the unit cell. We use angstrom as our unit of choice.

scf K_POINTS input

K_POINTS automatic
  9 9 1 1 1 1

Description

The K_POINTS card sets the k points at which QE calculates the wavefunctions and energies. The automatic option with input $n_1, n_2,n_3,s_1, s_2, s_3$ generates a Monkhorst-Pack grid of dimension $n_1\times n_2\times n_3$ with respective offsets $s_i=0\vee1$, corresponding to no offset or half a gridpoint offset.

The bands and nscf files are identical save for the calculation and K_POINTS cards:

nscf K_POINTS input

Description

The K_points card is list of points along a path in k space for bands and a grid for nscf, preceded by the number of points, in this case 64. Here we replace automatic with crystal which reads the k points in crystal coordinates (relative to the reciprocal lattice vectors). The fourth number is the relative weight of the k point.

Command line interface¶

The command line interface for QE is quite simple:

pw.x < prefix.scf.in > prefix.scf.out

where prefix is graphene in our case.

We can parallelize QE with the mpi library:

mpirun -n #n pw.x < prefix.scf.in > prefix.scf.out

here the -n flag can differ per distribution (sometimes it is -np). #n is the number of processes at our disposal. We can then choose how to divide these processes with additional flags after the pw.x executable. The main flag we used is the -nk flag which divides the calculation in #nk pools which each calculates a fraction of the k points, allowing for perfect parallelization. For higher #nk the RAM usage also increases so we can not automatically just set #nk to the amount of k points.

Bands output¶

The main output we directly care about is the bands output which leads a summary of the input arguments and the symmetries of the lattice. Since we define all K points explicitly these will not be used by QE. When we choose an automatic grid QE significantly decreases the points it needs to evaluate. Near the end of the file all k points and energies of the bands are printed.

bands output

          k = 0.0000 0.0000 0.0000 (  5649 PWs)   bands (ev):

   -23.8181 -11.9014  -7.2766  -7.2766  -1.3477  -0.9995  -0.7246  -0.5157
    -0.0219   0.2830   1.0299   1.4288   2.4049   2.9177   4.0662   4.1555
     4.1555   4.7234   6.0004   6.8461   7.2901   8.1773   8.4570   9.2604
    10.5975  11.9548  13.2678  14.9152  16.2216  18.1415

Description

The first k point evaluated by QE and all corresponding energies in eV.

In Fig. 15 we plot the QE bands together with the TB bands from (22). The model works quite well for the valence part of the QE bands. Ostensibly we require NNN hoppings and beyond to capture the conduction band component more accurately.

Projwfc input¶

Before we start the wannierization process we need to know which bands we want to wannierize and which atomic orbital contribute to these bands. QE offers a program called projwfc (project wavefunctions) which does exactly this.

projwfc input

&PROJWFC
  prefix = 'graphene'
  outdir = '../out'
  DeltaE = 0.05
/

Description

The PROJWFC card is the only card for projwfc source files. The prefix and outdir should match the project. DeltaE is the step size for the dos calculation of the projected orbitals, which we do not need but must be set for the program to work.

Projwfc output¶

PROJWFC output: states

     state #   1: atom   1 (C  ), wfc  1 (l=0 m= 1)
     state #   2: atom   1 (C  ), wfc  2 (l=1 m= 1)
     state #   3: atom   1 (C  ), wfc  2 (l=1 m= 2)
     state #   4: atom   1 (C  ), wfc  2 (l=1 m= 3)
     state #   5: atom   2 (C  ), wfc  1 (l=0 m= 1)
     state #   6: atom   2 (C  ), wfc  2 (l=1 m= 1)
     state #   7: atom   2 (C  ), wfc  2 (l=1 m= 2)
     state #   8: atom   2 (C  ), wfc  2 (l=1 m= 3)

Description

The projwfc output file contains a list of the projected atomic orbitals. For graphene we have two sites: atom 1 and atom 2, two atomic orbital types: s-orbitals with azimuthal quantum number l=0 (wfc 1) and p-orbitals with azimuthal quantum number l=1 (wfc 2). For l=1 there are three different possibilities for the magnetic quantum number $m_l=-1, 0, 1$. For a spinless bandstructure this quantum number is reformulated to m=1 for the p$_z$ orbital, m=2 for the p$_x$ orbital and m=3 for the p$_y$ orbital. In total we have 8 orbitals.

The second part of the projwfc output file list the contribution of each state to each band at each k points.

PROJWFC output: states

 k =   0.0000000000  0.0000000000  0.0000000000
==== e(   1) =   -23.81812 eV ==== 
     psi = 0.492*[#   1]+0.492*[#   5]+
    |psi|^2 = 0.984
==== e(   2) =   -11.90144 eV ==== 
     psi = 0.473*[#   2]+0.473*[#   6]+
    |psi|^2 = 0.946
==== e(   3) =    -7.27665 eV ==== 
     psi = 0.249*[#   3]+0.249*[#   4]+0.249*[#   7]+0.249*[#   8]+
    |psi|^2 = 0.997
==== e(   4) =    -7.27665 eV ==== 
     psi = 0.249*[#   3]+0.249*[#   4]+0.249*[#   7]+0.249*[#   8]+
    |psi|^2 = 0.997
==== e(   5) =    -1.34766 eV ==== 
     psi = 0.002*[#   1]+0.002*[#   5]+
    |psi|^2 = 0.005
==== e(   6) =    -0.99950 eV ==== 
     psi = 0.004*[#   2]+0.004*[#   6]+
    |psi|^2 = 0.007

Description

The contribution of the eight states to the first 6 bands at the first k point. If we write the wavefunction as $\Psi = \sum\alpha_i\psi_i$ then the number before each state corresponds to $|\alpha_i|^2$. We see that for bands 5 and 6 |psi|^2 is not even close to 1. This tells us that this band can not be described with these 8 states. The leftover plus at the end of the line is just a for loop artifact which can be ignored.

In tandem with the projected orbitals we wrote a small script to untangle the bands so we can separate the future Wannier bands from the other conduction bands.

Wannier90¶

Next up is the wannierization of the two bands forming the Dirac cone. FromFig. 16 we learn that only the p$_z$ orbitals contribute and that these bands extend far into the valence bands. To untangle these bands we need to make sure that takes all these bands into account. That is why we have set nbnb quite high for the nscf calculation so that QE produces the whole band. It is essential that the nscf calculation is run after a bands calculation using the same outdir for both as they overwrite each other’s wavefunctions.

W90 input¶

W90 general input

guiding_centres = true
write_hr = true
wannier_plot = true
bands_plot = true
num_bands = 30
num_wann = 2
num_iter = 2000
wannier_plot_supercell = 3, 3, 1

Description

guiding_centres tells W90 to use the projection (orbital) centres as the guiding centres for the Wannier orbitals
write_hr tells W90 to write the hopping energies to a file,band_plot does the same for the Wannier bands and wannier_plot for the Wannier orbitals. We use wannier_plot_supercell = 3, 3, 1 to prevent Wannier orbitals to be cut off at the edges of the unit cell (extends the unit cell in the xy plane).
num_bands describes the number of bands coming from the QE nscf calculation and num_wann the number of Wannier bands to project to.
num_iter is the maximum number of iterations W90 is allowed to execute while minimizing the spread of the Wannier orbitals. If the spread is not satisfactory after the given number of iteration, add restart=wannierise to the input and W90 will continue where it left off.

W90 disentanglement input

dis_win_min = -12
dis_win_max = 8
dis_froz_min = -7
dis_froz_max = -1.5
dis_num_iter = 500

Description

These arguments set the parameters with which to disentangle the #num_bands to extract #num_wann bands.

The choice of dis_froz_min and dis_froz_max is the most important challenge of any W90 calculation. Together they set a window in which only the Wannier bands are present. The window should then be chosen as large as possible for better convergence. In Fig. 16 we clearly see that between energies -7 eV and -1.5 eV only the Dirac bands are present, an easy choice for the disentanglement window.
The dis_min and dis_max arguments also form a window which should contain all Wannier orbitals and be as small as possible. Fig. 16 indicates that a good choice is the window from -12 eV to 8 eV.
dis_num_iter is de maximum amount of disentanglement iterations to perform.

W90 lattice input

begin unit_cell_cart
  2.4595121   0.0000000   0.0000000
  -1.229756   2.1299999   0.0000000
  0.0000000   0.0000000   32.000000
end unit_cell_cart

Begin atoms_cart
ang
  C   0.0000000   0.0000000   0.0000000
  C   0.0000000   1.4200000   0.0000000
End atoms_cart

Begin projections
  C:  pz
End projections

Description

These arguments set the lattice and projections.

unit_cell_cart lists the three lattice vectors in Angstrom.
atoms_cart lists all atoms and their location in the unit cell.
projections sets the initial guess for the Wannier orbitals as some mix of atomic orbitals. Fig. 16 tells us that a good guess for the Dirac bands are the p$_z$ orbitals. W90 automatically adds this projection for every atom of the same type. It is also possible to add atomic orbitals to specific locations if we want more control of the exact placement. It is essential that the total amount of projections is identical to num_wann.

W90 band_plot input

Begin Kpoint_Path
  Γ 0.0000000   0.0000000   0.0000000   M 0.5000000   0.0000000   0.0000000
  M 0.5000000   0.0000000   0.0000000   K 0.3333333   0.3333333   0.0000000
  K 0.3333333   0.3333333   0.0000000   Γ 0.0000000   0.0000000   0.0000000
End Kpoint_Path

Description

The Kpoint_Path is the path in k space W90 follows for the band_plot and has no direct link to the QE calculation.

W90 kpoints input

Description

These k points have to be identical to the k points fo the nscf calculation which is why it is essential to manually define them in the QE input file. mp_grid sets the dimensions of the Monkhorst-Pack grid which has to match the total number of k points.

PW2WAN input¶

The PW2WAN program is part of the QE code which exports the QE wavefunctions for postprocessing by W90.

PW2WAN input

&inputpp
   outdir = './out'
   prefix = 'graphene'
   seedname = 'graphene'
   write_unk = .true.
   write_mmn = .true.
   write_amn = .true.
/

Description

The outdir and prefix should match the outdir and prefix of the QE calculation whereas the seedname should match the name if the W90 file (minus the extension .win).
write_mmn and write_amn are crucial inputs for W90 whereas write_unk is only necessary for plot_wannier. For a large amount of bands the unk files can become massive (couple of GB a piece times the number of k points) so we advise to only add this parameter for small lattices.

Command line interface¶

Before we can run PW2WAN we need to preprocess the W90 input file. So we run the following series of commands:

wannier90.x -pp seedname.win
pw2wannier90.x < prefix.pw2wan.in > prefix.pw2wan.out
wannier90.x seedname.win

W90 offers an mpi version for parellization which has to be compiled seperately. We can parallelize pw2wannier90.x in the same way as pw.x.

W90 output¶

W90 disentanglement output

 +---------------------------------------------------------------------+<-- DIS
 |  Iter     Omega_I(i-1)      Omega_I(i)      Delta (frac.)    Time   |<-- DIS
 +---------------------------------------------------------------------+<-- DIS
       1       1.98367895       1.96460709       9.708E-03      0.00    <-- DIS
       2       1.97124443       1.95760453       6.968E-03      0.00    <-- DIS
       3       1.96331159       1.95227135       5.655E-03      0.00    <-- DIS

\[\vdots\]

     258       1.91625911       1.91625911       9.655E-11      0.73    <-- DIS
     259       1.91625911       1.91625911       9.135E-11      0.73    <-- DIS
     260       1.91625911       1.91625911       8.644E-11      0.73    <-- DIS

             <<<      Delta < 1.000E-10  over  3 iterations     >>>
             <<< Disentanglement convergence criteria satisfied >>>

        Final Omega_I     1.91625911 (Ang^2)

Description

Omega_I is the gauge invariant part of the spread which minimizes for smooth Wannier bands. Delta is the fractional change of Omega_I:

$$\Delta = \frac{|\Omega_I(i)-\Omega_I(i-1)|}{\Omega_I(i)}$$.

When this value drops below the convergence threshold dis_conv_tol the disentanglement process stops.

W90 wannierization output

 *------------------------------- WANNIERISE ---------------------------------*
 +--------------------------------------------------------------------+<-- CONV
 | Iter  Delta Spread     RMS Gradient      Spread (Ang^2)      Time  |<-- CONV
 +--------------------------------------------------------------------+<-- CONV

 *----------------------------------------------------------------------------*
 Initial State
  WF centre and spread    1  ( -0.000000,  0.000000,  0.000000 )     0.96825063
  WF centre and spread    2  ( -0.000000,  1.420000,  0.000000 )     0.96825057
  Sum of centres and spreads ( -0.000000,  1.420000,  0.000000 )     1.93650120

      0     0.194E+01     0.0000000000        1.9365012594       0.73  <-- CONV
        O_D=      0.0133595 O_OD=      0.0068826 O_TOT=      1.9365013 <-- SPRD
 *----------------------------------------------------------------------------*

\[\vdots\]

 *----------------------------------------------------------------------------*
 Cycle:   2000
  WF centre and spread    1  (  0.000000,  0.000000,  0.000000 )     0.96817487
  WF centre and spread    2  ( -0.000000,  1.420000,  0.000000 )     0.96817481
  Sum of centres and spreads (  0.000000,  1.420000,  0.000000 )     1.93634968

   2000     0.000E+00     0.0000000100        1.9363497402       1.85  <-- CONV
        O_D=      0.0132080 O_OD=      0.0068826 O_TOT=      1.9363497 <-- SPRD
 Delta: O_D=  0.1734723E-16 O_OD= -0.1734723E-17 O_TOT=  0.0000000E+00 <-- DLTA
 *----------------------------------------------------------------------------*

 Writing checkpoint file graphene.chk... done

 Final State
  WF centre and spread    1  (  0.000000,  0.000000,  0.000000 )     0.96817487
  WF centre and spread    2  ( -0.000000,  1.420000,  0.000000 )     0.96817481
  Sum of centres and spreads (  0.000000,  1.420000,  0.000000 )     1.93634968

         Spreads (Ang^2)       Omega I      =     1.916259106
        ================       Omega D      =     0.013208013
                               Omega OD     =     0.006882621
    Final Spread (Ang^2)       Omega Total  =     1.936349740

Description

In the wannierization process the gauge dependent spread (Omega_D+Omega_OD, O(D) for (off)-diagonal) is minimized with a steepest descend algorithm ( RMS gradient is the root mean square of the gradient). Delta does not drop below the convergence threshold so W90 completes all num_iter cycles. Our projections work quite well for this system so not much optimization is needed.

The NN hopping energy between two Wannier orbitals is approximately 2.9 eV, slightly more than the 2.7 eV of [CDAnjouG+11], which aims to describe the Dirac point more accurately.

The Wannier orbitals of Fig. 18 closely resemble the p$_z$ projections, as expected from the projwfc calculation in Fig. 16.

Wannierization WSe$_2$

hidde-dijkstra.github.io

Graphene¶

Basic theory¶

Quantum Espresso¶

SCF/NSCF/Bands input¶

Command line interface¶

Bands output¶

Projwfc input¶

Projwfc output¶

Wannier90¶

W90 input¶

PW2WAN input¶

Command line interface¶

W90 output¶